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INTRODUCTION: 


Effective  computational  tools  are  developed  for  improved  detection,  diagnosis,  and 
monitoring  of  breast  cancer.  A  tactile  mapping  device  (TMD)  is  tested  for  the  e^ly 
detection  of  breast  cancer  through  more  objective  and  quantitative  breast  palpation. 
Following  the  detection  of  breast  cancer,  diagnosis  of  the  breast  is  the  next  step.  We 
introduce  a  hybrid  source  decomposition  algorithm,  which  allows  for  a  computational 
characterization  of  tumor  microvascular  heterogeneity  in  both  spatial  and  temporal 
domains.  The  goal  is  to  reveal  temporal-spatial  patterns  for  the  visualization  and 
quantification  of  tumor-induced  angiogenesis  and  response  to  therapy  using  dynamic 
contrast-enhanced  magnetic  resonance  imaging  (DCE-MRI).  An  image-based  change 
detection  approach  is  proposed  to  assess  tumor’s  response  to  therapy.  Based  on 
promising  experimental  results,  we  anticipate  that  fimctional  imaging  based 
computational  characterization  of  tumor  heterogeneity  and  response  will  be  useful  in  a 
wide  variety  of  medical  imaging  studies. 

Our  considerable  efforts  in  the  past  year  are: 

•  To  construct  training  and  testing  databases  of  tactile  mapping  images  and 
associated  data. 

•  To  construct  MR  modeling  of  tumor  to  estimate  and  tract  changes  of  the  tumor 
across  time. 

•  To  evaluate  the  performance  of  the  TMD+NN  in  diagnosis  and  monitoring  in 
terms  of  sensitivity  to  changes  and  reproducibility  to  measurements. 
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BODY: 


In  the  first  year  of  this  award,  we  have  fully  demonstrated  the  feasibility  of  the  tactile 
mapping  device  (TMD)  for  improving  breast  cancer  examination  technique  in  diagnosis, 
documentation,  and  training.  In  particular,  the  results  have  shown  that  new  tactile 
mapping  technology  can  quantitatively  measure  the  location  and  applied  forces  in  breast 
palpation,  and  the  tactile  features  of  detected  breast  lumps.  The  prototype  interactive 
training  program  can  track  finger  motions  and  applied  forces  during  breast  palpation  in 
which  on-line  feedback  can  help  the  training  to  better  understand  the  search  strategy  and 
adjust  applied  force  level  to  increase  the  sensitivity. 

In  the  second  year  of  this  award,  we  have  integrated  the  tactile  sensing  technology  and 
the  vision-based  neural  network  for  investigations  of  soft  tissue  interactions  with  the 
tactile/force  sensor.  With  the  proven  power  of  nonlinear  signal  processing  both 
convolution  neural  network  (CNN)  and  multi-layer  perceptron  (MLP)  can  be  sued  to 
characterized  the  hard  inclusion  (breast  cancer)  through  neural  network  learning 
capabilities,  instead  of  using  a  simplified  complex  biomechanics  model  with  many 
heuristic  assumptions.  The  tactile  mapping  systems  using  the  neural  networks  and  tactile 
sensing  array  can  extract  invariant  parameters  associated  with  the  lesions  (i.e.,  the  size 
and  depth  of  the  lesion)  can  be  estimated  more  accurately  than  those  by  the  conventional 
approaches. 

Following  the  detection  of  breast  cancer,  diagnosis  of  the  breast  is  the  next  step. 
Recently,  there  has  been  a  need  to  stimulate  the  development  of  novel  imaging 
technologies  that  exploit  our  current  knowledge  of  the  genetic  and  functional  bases  of 
environmentally  induced  diseases  and  cancers.  Such  functional  imaging  capabilities  will 
allow  for  the  visualization  and  explanation  of  important  disease-causing  physiological 
and  functional  processes  in  the  living  tissue.  Therefore,  the  potential  of  medical  imaging 
to  improve  cancer  treatment  extends  well  beyond  using  imaging  information  to  help 
select  effective  preventatives  or  treatments. 

Dynamic  contrast-enhanced  magnetic  resonance  imaging  (DCE-MRJ)  has  been 
developed  to  provide  additional  functional  information  about  tumor  [1].  The  technique 
has  emerged  as  an  effective  tool  to  access  tumor  vascular  characteristics,  which  can  be 
used  for  measuring  tissue  perfusion,  microvessel  permeability,  and  vascular  volume.  The 
imaging  technique  has  provided  the  information  needed  for  diagnosis  and  treatment 
based  on  tumor  location,  size,  and  spread.  It  is  also  known  that  the  kinetic  characteristic 
changes  following  treatment  have  correlated  with  histopathological  outcome  (e.g.,  tumor 
blood  volume,  vascular  permeability,  and  tumor  perfusion)  and  patient  survival.  However 
pathology  reveal  distinct  heterogeneity  within  and  among  tumors,  that  suggesting  the 
need  for  individualize  diagnosis  and  treatment  to  the  unique  characteristics  of  each 
specific  case.  Thus,  dynamic  measurements,  in  which  the  upt^e  and  washout  of  contrast 
in  tissues  is  monitored  with  time,  can  assist  in  the  diagnosis  of  breast  tumors  and  can 
provide  information  on  vascular  permeability  and  perfusion.  The  techniques  are  being 
applied  to  assess/monitor  the  response  to  treatment,  which  can  be  used  to  characterize 
microvasculature  providing  information  about  tumor  microvessel  structure  and  function. 
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However,  widespread  success  of  DCE-MRI  may  be  limited  by  the  need  for  the  further 
development  of  technology,  particularly  due  to  the  lacking  of  quantitative  and 
computational  data  analysis  tools  included  by  the  instruments  [1].  Therefore,  an  effective 
computational  data  analysis  technique  for  image-based  lesion  characterization  of  the 
vascularlity  patterns  (fasWslow-flow  kinetic)  in  functional  imaging  will  provide  richer 
diagnostic  information,  thus  improve  the  sensitivity  (e.g.,  signal-to-noise  ratio)  and 
specificity  (e.g.,  heterogeneity)  for  early  disease  detection/diagnosis  and  monitoring  of 
therapeutic  response.  The  techniques  are  mainly  used  for  the  prognosis  in  breast  cancer, 
and  it  also  shows  more  promising  for  assessing  the  response  to  therapy.  Subsequently, 
functional  imaging  will  play  an  important  role  in  the  early  detection,  diagnosis,  and 
treatment  of  diseases  [2]. 

In  this  research,  we  introduced  a  hybrid  decomposition  algorithm,  which  allows  for  a 
computed  simultaneous  imaging  of  multiple  biomarkers.  The  method  is  based  on  a 
combination  of  time-activity  curve  clustering,  pixel  subset  selection,  and  independent 
component  analysis.  We  demonstrate  the  principle  of  the  approach  on  an  image  data  set, 
and  we  then  apply  the  method  to  the  tumor  vascular  characterization  using  DCE-MRI. 

Major  portions  of  the  work  were  reported  in  our  manuscripts  submitted  to  the 
SPIE’s  medical  imaging  conference  and  NNSP  2003  conference.  [Attached] 

KEY  RESEARCH  ACCOMPLISHMENTS: 

Our  key  goal  of  Task  11  is  to  develop  a  neural  network  based  intelligence  system  that 
estimates  and  track  the  changes  (in  size  and  depth)  of  tumors  over  time  after  diagnosis 
and  during  treatment.  Furthermore,  the  MR  modeling  of  tumor  to  estimate  and  track 
changes  of  the  tumor  across  time  in  diagnosis  and  treatment  is  presented.  Our 
accomplishments  are  the  following: 

•  We  have  developed  neural  network  training  algorithms  and  computer 
interpretation  codes,  convolution  neural  network  (CNN)  based  breast  tumor 
characterization  and  parameter  estimation. 

•  We  have  constructed  the  databases  of  tactile  mapping  images  and  the  associated 
data  for  the  training  and  testing  of  our  neural  network  algorithms. 

•  We  have  developed  a  hybrid  source  decomposition  algorithm,  which  allows  for  a 
computational  characterization  of  tumor  microvascular  heterogeneity  in  both 
spatial  and  temporal  domains.  The  goal  is  to  reveal  temporal-spatial  patterns  for 
the  visualization  and  quantification  of  tumor-induced  angiogenesis  and  response 
to  therapy  using  dynamic  contrast-enhanced  magnetic  resonance  imaging  (DCE- 
MRI). 


3 


REPORTABLE  OUTCOMES: 


•  R  Srikanrhana  .  Functional  Imaging  and  Analysis  of  Tumor  Heterogeneity  by 
Cluster  and  Independent  Component  Analysis,  Doctoral  dissertation,  The  Catholic 
University  of  America,  2003 . 

•  Y.  Wang,  J.  Zhang,  R.  Srikanchana.  J.  Xuan,  Z.  Wang,  Z.  Szabo,  Z.  Bhujwalla,  P . 
Choyke,  K.  Li,  “Computed  Simultaneous  Imaging  of  Multiple  Biomarkers,”  Proc. 
IEEE  Neural  Networks  for  Signal  Processing,  2003. 

•  J.  Zhang,  R.  Srikanchana.  J.  Xuan,  K.  Li,  Y.  Wang  “Partially-Independent 
Component  Analysis  for  Molecular  Imaging,”  SPIE’s  Inti.  Symp.  Medical 
Imaging,  vol.  5032,  San  Diego,  CA,  February  2003 

•  Presented  the  paper  “Partially-Independent  Component  Analysis  for  Molecular 
Imaging,”  SPIE's  Inti  Symp  Medical  Imaging,  San  Diego,  CA,  Feb.  2003. 

CONCLUSIONS: 

In  this  research,  effective  computational  tools  are  developed  for  improved  detection, 
diagnosis,  and  monitoring  of  breast  cancer  and  response  to  therapy.  We  have 
demonstrated  the  principle  of  the  approaches  on  both  simulated  and  real  data  sets.  Based 
on  promising  experimental  results,  we  anticipate  that  functional  imaging  based 
computational  characterization  of  tumor  heterogeneity  and  response  will  be  useful  in  a 
wide  variety  of  medical  imaging  studies. 

The  main  contribution  of  this  research  have  shown  that  the  TMD  with  associated  new 
tools  will  improve  physical  breast  examination  in  the  ability  to  detect  small  tumors  in 
breast  palpation  and  increase  the  accuracy  of  early  detection,  leading  to  improved 
diagnosis  and  a  reduction  in  breast  cancer  mortality.  We  have  adapted  new  tactile 
mapping  technology  to  the  needs  of  improving  physical  breast  examination  and  gather 
preliminary  clinical  data.  We  have  conducted  a  study  to  advance  fundamental 
understanding  of  palpation  and  solve  these  practical  problems  through  the  creation  of  a 
new  TMD.  This  device  measured  three  key  variables  during  palpation:  the  examiners 
search  patterns,  the  applied  forces,  and  the  small-scale  pressure  variations  at  the  skin  due 
to  lumps.  A  preliminary  TMD  prototype  demonstrated  the  feasibility  and  effectiveness. 
We  have  heavily  focused  on  the  DCE-MRI  based  angiogenesis  imaging  to  elucidate  two 
main  challenging  aspects  (tumor  heterogeneity  and  composite  signals)  in  data  analysis  of 
functional  imaging.  Complementary  to  various  existing  methods  (e.g.,  compartment 
modeling,  factor  analysis),  we  introduced  a  hybrid  source  decomposition  algorithm, 
which  allows  for  a  computational  characterization  of  tumor  microvascular  heterogeneity 
in  both  spatial  and  temporal  domains.  The  method  is  based  on  a  combination  of  time- 
activity  curve  clustering,  pixel  subset  selection,  and  independent  component  analysis. 
The  goal  is  to  reveal  temporal-spatial  patterns  for  the  visualization  and  quantification  of 
tumor-induced  angiogenesis  and  response  to  therapy  using  DCE-MRI  of  breast  cancer. 
As  a  result,  spatial  distribution  of  tumor  blood  volume,  vascular  permeability,  and  tumor 
perfusion,  as  well  as  their  TACs  are  simultaneously  estimated,  which  closely  resemble 
the  expected  characteristics  of  the  tumor  heterogeneity.  The  method  can  be  used  to 
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evaluate  tumor  angiogenesis  inhabitation  for  cancer  therapy  and  monitor  the  kinetic 
characteristic  response  to  the  anti-angiogenic  following  the  treatment.  The  method 
consists  of  two  major  steps:  (1)  a  multivariate  cluster  analysis  for  the  initialization  of  the 
factor  image  decomposition;  (2)  a  factor  image  decomposition  by  a  partially-independent 
component  analysis 

We  wish  to  suggest  that  our  preliminary  studies  indicated  a  promising  utility  of  hybrid 
blind  source  separation  techniques  in  computed  simultaneous  imaging  of  multiple 
biomarkers.  Although  the  optimality  of  this  method  may  be  data-dependent,  we  would 
expect  it  to  be  an  effective  tool  applicable  to  various  modalities. 
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Abstract.  Functional- molecular  imaging  techniques  promise  pow¬ 
erful  tools  for  the  visualization  and  elucidation  of  important  disease- 
causing  physiologic- molecular  processes  in  living  tissue.  Most  ap¬ 
plications  aim  to  find  temporal— spatial  patterns  assocaited  with 
different  disease  stages.  When  multiple  agents  are  used,  imagery 
signals  often  represent  a  composite  of  more  than  one  distinct  source 
due  to  functional-molecular  biomarker  heterogeneity,  independent 
of  spatial  resolution.  We  therefore  introduce  a  hybrid  decomposi¬ 
tion  algorithm  which  allows  for  a  computed  simultaneous  imaging 
of  multiple  biomarkers.  The  method  is  based  on  ^  combination  of 
time-activity  curve  clustering,  pixel  subset  selection,  and  indepen¬ 
dent  component  analysis.  We  demonstrate  the  principle  of  the  ap¬ 
proach  on  an  image  data  set,  and  we  then  apply  the  method  to  the 
tumor  vascular  characterization  using  dynamic  contrast- enhanced 
magnetic  resonance  imaging  and  brain  neuro- transporter  imaging 
using  dynamic  positron  emission  tomography. 


INTRODUCTION 

With  rapid  advances  recently  made  in  developing  molecular /functional-targeted 
contrast  agents,  ligands  and  imaging  probes,  new  molecular  or  functional 
imaging  techniques  promise  powerful  tools  for  the  visualization  and  eluci- 
.dation  of  important  disease-causing  physiologic  and  molecular  processes  in 
living  tissue  [1,  2].  For  example,  dynamic  contrast-enhanced  magnetic  reso¬ 
nance  imaging  (DCEkMRI)  utilizes  various  molecular  weight  contrast  agents 
to  access,  non-invasively,  tumor  microvascular  characteristics.  The  extravas- 


cular  retention  of  intravenous  contrast  medium  correlates  to  accumulation  at 
sites  of  concentrated  angiogenesis  mediating  molecules  or  microvessels  (per¬ 
meability).  Kinetic  (perfusion)  changes  following  treatment  have  correlated 
with  histopathological  outcome  and  patient  survival,  and  shall  potentially  be 
able  to  assess  or  predict  the  response  to  treatment  particularly  when  using 
anti- angiogenic  drugs  [2].  On  the  other  hand,  positron  emission  tomography 
(PET)  utihzes  molecular  probes  (e.g,  hgands  for  receptors  or  substrates  for 
intracellular  enzymes)  labeled  with  positron-emitting  radioisotopes.  Tracer 
quantities  yield  a  tomographic  image  after  their  retention,  as  a  consequence  of 
either  specific  ligand-receptor  binding  or  conversion  of  substrate  to  “trapped” 
metabolic  product (s)  [1].  Fig.  1  shows  the  breast  tumor  images  obtained  by 
DCE-MRI. 


Figure  1:  Images  of  advanced  breast  tumor  obtained  by  DCE-MRI,  where  tumor- 
induced  microvascular  permeabilities,  associated  with  different  perfusion  rate,  re¬ 
veal  interesting  spatial  heterogeneous  patterns. 

As  a  common  feature  in  functional  or  molecular  imaging,  pixels  repre¬ 
sent  a  composite  of  more  than  one  distinct  biomarker,  independent  of  the 
spatial  resolution.  For  example,  DCE-MRI  reveals  heterogeneous  mixture 
of  tumor  microvessels  associated  with  different  perfusion  rate  (e.g.,  fast  vs. 
slow  flow),  and  PET  imaging  of  the  neuro- transporters  in  the  brain  indicates 
that  although  highly  sensitive  to  its  target,  the  nonspecific  binding  of  the 
radioligand  is  significant.  As  a  result,  the  overlap  of  multiple  biomarkers 
can  severely  decrease  the  sensitivity  and  specificity  for  the  measurement  of 
functional  or  molecular  signatures  associated  with  different  disease  processes. 

Two  most  popular  methods  that  aim  to  extract  the  images  of  individual 
biomarkers  are  compartment  modeling  and  factor  analysis  [3j.  Fig.  2  shows 
macroscale  kinetics  of  probe-biomarker  interactions.  Compartment  model¬ 
ing  can  find  time-activity  curves  (TACs)  which  are  functionally  meaningful, 
but  requires  pre-acquired  input  function  and  does  not  generally  provide  any 
spatial  information  about  different  tissue  kinetics  [3].  Factor  analysis  can 
find  both  factor  images  and  TACs  including  the  input  function  blindly,  but 
may  result  in  the  factor  images  or  TACs  that  are  not  unique  or  functionally 
meaningful  [4].  Several  blind  compartment  modeling  methods  were  recently 
proposed,  where  the  major  concern  is  the  parameter  identifiabihty  ■  [5] . 


Figure  2:  Dlustration  of  two-tissue  compartmentl  model  and  time-activity  curves. 

In  this  paper,  we  introduce  a  hybrid  decomposition  algorithm  which  allows 
for  a  computed  simultaneous  imaging  of  multiple  biomarkers.  The  method  is 
based  on  a  combination  of  time- activity  curve  clustering,  pixel  subset  selec¬ 
tion,  and  independent  component  analysis.  We  demonstrate  the  principle  of 
the  approach  on  an  image  data  set,  and  we  then  apply  the  method  to  the  tu¬ 
mor  vascular  characterization  using  DCE-MRI  and  brain  neuro- transporter 
imaging  using  dynamic  PET. 


THEORY  AND  METHOD 


Factored  Compartment  Model 


We  first  introduce  a  simple  form  of  linear  factored  compartment  model  and 
discuss  its  application  to  simultaneous  biomarker  imaging.  On  the  basis  of 
Fig.  2,  tracer  characterization  within  a  region  of  interest  (ROI)  leads  to  a 
set  of  first  order  differential  equations: 


Cf{t)  =  kifCj,(t)  -  k2fCf{t)  ,  . 

Csit)  =  kuCplt)  -  k2sCsit)  , 

c^(t)  =  Cf(t)+Csit)  +  Cpit) 


fciyCp(t)  0  6-^2/* 
kisCp(t)  ® 


where  Cf{t)  and  Cs{t)  are  the  tissue  activity  in  the  fast  turnover  and  slow 
turnover  pools,  respectively,  at  time  t;  Cp(t)  is  the  tracer  concentration  in 
plasma  (i.e.,  the  input  function);  Cm(t)  is  the  measured  total  tissue  activity; 
kjf  and  kis  are  the  unidirectional  transport  constants  from  plasma  to  tissue 
(permeability  in  ml/min/g:  spatially-var5dng);  k2f  and  k2s  the  rate  con¬ 
stants  for  efflux  (perfusion  in  /min:  spatially-invariant) ,  and  0  denotes  the 
mathematical  convolution.  Note  that  c/(t),  Cs{t)j  Cp(t)^  and  CTn{t)  are  called 
the  time-activity  curves  (TACs). 

Linear  system  theory  suggests  a  simple  method  to  convert  temporal  kinet¬ 
ics  to  spatial  information  [6].  Let  kif{i)  and  kis{i)  be  the  local  permeability 
parameters  associated  with  pixels  i  -  i.e.,  the  permeability  of  fast 

and  slow”  turnover  (perfusion)  regions  in  the  pixel,  respectively,  and  Vp{i)  be 
the  plasma  volume  in  tissue.  The  measured  pixel  TAG  Cm[i^t)  is  [2] 


=  kif(i)Cp{t)  ®  +  ku{i)cp(t)  ®  e  +  Vp(i)cp{t)  (2) 


that  leads  to  a  factored  compartment  model  [4] 

1  r  Cp(ti)  ®  Cpih)  ®  Cp(ti)  ' 

CmiiM)  Cp(t2)  ®  Cp(f2)  ®  Cp(t2) 

=  ■  ■  ■  ■  kxs{i) 

L 

^n)  ^pi^n)  ^2f^n  Cp(tn)  ®  e  "  Cp(tn) 

(3) 

where  A=  [  Cp(f)  ®  c^(i)  ®  Cp(i)  ]  is  called  mixing  matrix 

consisting  of  the  TACs  associated  with  underlying  factor  images  whose  least 
square  estimates  are 

ki.ii)  =  (a'^a)  A^  •  .  (4) 

vp(i)  J  : 

C7n(^)^n) 

Below  we  propose  three  different  approaches  that  aim  to  reconstruct  the 
factor  images  for  the  cases  of  multi-region  insignificantly-overlapped,  single- 
region  significantly-overlapped,  and  multi-region  significantly-overlapped,  het¬ 
erogeneities. 

Multivariate  Cluster  Analysis 

Biomarker  heterogeneity  in  space  motivates  a  natural  consideration  of  mul¬ 
tivariate  cluster  analysis  (i.e.,  pixel  TAG  classification)  for  factor  image  de¬ 
composition.  Since  vascular  space  Cp(t)  is  usually  very  high  in  the  early 
time-course  and  dropped  quickly  (e.g.,  delta  function),  our  initial  analysis 
focuses  on  the  tissue  perfusion  activities.  Ideally,  when  there  are  only  pure- 
volume  pixels,  i.e.,  [  kif(i)  kis{i)  has  only  one  nonzero  element,  there 
shall  be  a  one-to-one  mapping  between  pixel  TAG  Cm(^^)  ^.nd  factor  TAG 
Cf{t)  or  Cs(t)  except  for  an  arbitrary  scaling  and  noise  eflfect. 

Since  the  shape  rather  than  the  magnitude  of  the  pixel  TAG  is  of 

the  interest,  we  perform  “centering”  and  “normalization”  on  each  pixel  TAG 
Cm  (^5^)  over  time  t  to  a  constant  scale  with  mean  zero 

1  "" 

Cm  (^*5  i)  ””  “  )  /^<^7n  (i)  (^) 

where  <TcTn(0  standard  de\dation  of  Cm(i,t)  over  t.  We  have 

t) = Kimit) + kUi)c*At),  ktfii) + ki(i)  =  1  (6) 

where  each  of  can  take  on  any  value  in  the  range  (0, 1)  subject 

to  the  noise  effort.  Assume  that  the  noise  effect  is  approximately  gaussian, 


the  pixel  TAG  c^{i,  t)  can  be  represented  by  a  gaussian  random  vector  (z) 
with  mean  fxj  ^  and  covariance  matrix  There  has  been  considerable 

success  in  using  the  standard  finite  normal  mixture  (SFNM)  distribution  to 
adequately  model  clustered  data  set,  taking  a  sum  of  the  following  general 
form  [7] 

P(Cm)=’r/5«lM/,C/)+7r,p(c*  |ax,,C,)  ^  (7) 

where  TTf^g  is  the  corresponding  mixing  proportion,  with  0  <  <  1  and 

TTy  -f  TT^  =  1,  and  g  is  the  gaussian  kernel. 

The  maximum  likelihood  estimation  of  the  SFNM  model  can  be  performed 
by  the  expectation-maximization  (EM)  algorithm.  We  have  previously  de¬ 
veloped  a  Visual  Data  Analyzer  (VISDA)  to  perform  multivariate  cluster 
analysis  [7,  8].  The  procedure  using  VISDA  provides  “soft”  sphts  of  the 
pixel  TAG  vectors,  hence  allowing  the  pixel  TAGs  to  contribute  si¬ 

multaneously  to  the  multiple  factor  TAGs.  Specifically,  the  outcome  of  such 
a  multivariate  cluster  analysis  is  the  posterior  Bayes  probabilities  of  pixel  i 
associated  with  c/,s(t) 

that  leads  to 

N  N  '  ' 

and  Cs{t)  =  Y^ZisCmii,t).  (9) 

i=:l  i=l 

Having  determined  the  factor  TAGs  A=  [  Cf(t)  Cs{t)  ],  the  factor  im¬ 
ages  k=  [  kif{i)  kis(})  will  be  reconstructed  from  Cmiiyt)  according  to 
Eq.  4.  However,  pixels  called  partial-volume  pixels  represent  a  combination 
of  factor  TAGs  due  to  the  partial-overlap  of  the  factor  images  in  space.  The 
partial-volume  effect  can  be,  fortunately,  solved  by  an  appropriate  partial- 
volume  modeling  [9].  Let  Cfs(t)  be  the  TAG  of  the  partial- volume  pixels. 
Based  one  the  derivation  in  [9],  we  can  write  a  SFNM  model  that  includes 
partial-volume  effect  as 

+T^fs9  (cm  I5  (M/  +  Ms)  )  J  (C/  +  Cs)  +  -^  -  fJ.^)  {flf  -  Ms)  ) 

(10) 

where  0  <  <  1  and  tt^  -f  -f  =  1.  We  have  modified  the  VISDA 

algorithm  to  estimate  this  constrained  SFNM  model.  We  will  first  apply 
unconstrained  EM  algorithm  to  classify  pixel  TAGs.  into  three  clusters  and 
then,  identify  two  clusters  with  smallest  covariances  as  those  representing 
pure- volume  pixels,  and  lastly  re-estimate  the  model  using  a  constrained  EM 
algorithm.  In  M-step  specifically,  we  will  only  estimate  the  parameters  of 
pure- volume  clusters  using  ^i(/,s),  and  assign  the  parameter  values  for  the 
partial- volume  cluster.  Once  again,  bsised  on  newly  estimated  factor  TAGs, 
the  factor  images  are  accordingly  reconstructed. 


Independent  Component  Analysis 


For  single-region  significantly- overlapped  cases,  the  similarity  between  Eq.  3 
and  latent  variable  modeling  motivates  the  consideration  of  independent  com¬ 
ponent  analysis  (ICA)  approach  [10].  As  aforementioned,  the  factor  images 
are  not  observable,  and  nothing  is  known  about  the  properties  of  the  TACs. 
In  the  absence  of  this  information,  one  has  to  proceed  “blindly”  to  recover  the 
factor  images  from  their  TAC-modulated  activity  mixtures  [4].  ICA  method 
utilizes  independency  as  a  guiding  principle  and  performs  a  nongaussian  fac¬ 
tor  analysis  leading  to  a  unique  solution  [11].  More  precisely,  by  assuming 
that  the  hidden  components  are  statistically  independent  with  nongaussian 
distributions,  these  hidden  sources  can  be  found  by  ICA,  except  for  an  ar¬ 
bitrary  scaling  of  each  signal  component  and  permutation  of  indices.  ICA 
approach  exploits  primarily  temporal  diversity  in  that  the  dynamic  images 
taken  at  different  times  carry  different  mixtures  of  the  factor  images  [4]. 

Specifically,  let  demixing  matrix  W  be  an  estimate  that  gives  a  good  ap¬ 
proximation  of  the  inverse  of  A,  then  according  to  the  Central  Limit  Theorem 
that  states  “a  sum  of  independent  random  variables  usually  has  a  distribu¬ 
tion  that  is  closer  to  gaussian  than  any  of  the  original  random  variables,”  the 
recovered  factor  image 

y(j)=Wc„(0  =  (WA)k(i)  (11) 

is  usually  more  gaussian  than  any  of  kj(i)  and  becomes  least  gaussian  when 
it  in  fact  equals  one  of  the  kj(i),  i.e.,  each  row  or  column  of  WA  has  only  one 
nonzero  element  where  W=PDA‘“^  with  P  and  D  being  permutation  and 
scaling  matrices.  Thus,  maximizing  the  nongaussianity  of  the  output  signal 
y(i)  produces  the  independent  components,  i.e.,  true  k(i). 

There  are  several  ways  for  estimating  the  model  of  ICA  including  maxi¬ 
mization  of  nongaussianity,  maximum  likelihood,  and  minimizing  joint  mu¬ 
tual  information.  Most  estimation  principles  and  objective  ftmctions  are 
equivalent,  at  least  in  theory.  Two  popular  and  pubhcly  available  software 
codes  are  FastICA  and  RunICA  algorithms.  A  thorough  discussion  on  ICA 
theory  can  be  found  in  a  recent  textbook  [10]. 

Partially-Independent  Component  Analysis 

For  multi-region  significantly- overlapped  cases,  we  have  foimd  that  direct  ap¬ 
plication  of  ICA  to  factor  image  decomposition  using  aU  the  pixels,  however, 
often  leads  to  an  unsatisfactory  recovery  of  k(z)  [4].  This  shall  not  be  a  sur¬ 
prise  since  a  more  reasonable  assumption,  suggested  by  underlying  anatomy, 
is  that  the  functional  regions  are  piecevdse  homogeneous.  This  seenaingly 
global  viewpoint  turns  out  to  have  important  consequences,  since  it  implies 
that  the  clustered  joint  distribution  of  the  factor  images  exhibits  only  lo-  . 
cal  independency  between  the  functional  regions  and  hence,  does  not  assure 
the  factor  images  to  be  globally  independent  [11].  Thus,  we  shall  expect  to 
achieve  a  better  decomposition  using  a  subset  of  pixels  that  corresponds  to 


a  homogeneously-overlapped  single-region  and  would  be  approximately  inde¬ 
pendent. 


Figure  3:  Simulated  image  sequence  consisting  of  ”tire"  and  "coin”  objects.  The 
mixing  matrix  resembles  the  typical  shapes  of  fast  and  slow  flows  in  DCE-MRI. 

We  therefore  developed  a  partiaUy-independent  component  analysis  tech¬ 
nique  [4].  Rather  than  using  all  the  pixels  that  give  rise  to  a  large  separation 
error,  we  attempt  to  identify  a  homogeneously-overlapped  single-region  and 
over  which  to  estimate  the  demixing  matrix  W  and  subsequently  factor  im¬ 
ages  k(i).  Once  again,  the  region  identification  can  be,  fortunately,  solved 
by  performing  multivariate  pixel-TAC  clustering  using  VISDA  software  [7,  8]. 
Note  that  in  the  present  case,  all  the  pixels  are  homogeneously-mixed  partial- 
volume  pixels.  Once  all  the  homogeneously-overlapped  individual  regions  are 
identified,  we  shall  perform  ICA  over  each  of  the  individual  pixel  subsets  and 
then  take  an  averaged  outcome  of  the  factor  images  and  TACs.  The  selection 
of  homogeneous  ROIs  can  be  done  by  either  visual  inspection  or  using  model 
selection  criteria  [7]. 

EXPERIMENTAL  RESULT 

We  shall  now  illustrate  the  operation  of  our  method  when  applied  to  a  simple 
image  data  set,  and  then  present  results  from  the  study  of  DCE-MRI  and 
PET  data  sets.  We  first  consider  a  digital  image  set  consisting  of  ten  ob¬ 
servations  generated  from  the  mixtures  of  “tire”  and  “coin”  images,  see  Fig 
3.  Experiments  are  conducted  over  several  runs  of  the  ICA  and  PICA  algo¬ 
rithms,  each  with  a  different  random  mixing  matrix.  With  an  appropriate 
pixel  selection,  the  estimated  TACs  are  given  against  the  true  TACs  in  Fig.  3 
(bottom  row).  In  this  experiment,  the  basic  principle  and  capable  nature  of 
the  PICA  method  are  evident  as  the  PICA,  using  only  independent  portion 


of  the  observations,  successfully  separates  two  image  patterns  that  have  been 
insufficiently  recovered  by  conventional  ICA  [4]. 

The  second  experiment  reports  the  effectiveness  of  multivariate  pixel  TAG 
clustering  for  multi-region  insignificantly- overlapped  cases.  The  data  set  is 
obtained  by  DCE-MRI,  see  in  Fig.  1.  This  is  an  advanced  breast  tumor 
case  that  shows  significant  vascular  heterogeneity  where  the  biomarkers  are 
vascular  permeability  with  fast  and  slow  perfusions.  The  two-dimensional 
projection  of  the  pixel  TACs  is  given  in  Fig.  4  (left)  and  the  estimated  fac¬ 
tor  TACs  are  given  in  Fig.  4  (right).  We  noticed  that  the  reconstructed 
factor  images  by  multivariate  cluster  analysis  paralleled  the  expected  vas¬ 
cular  distributions  very  nicely  and  many  independent  trials  reached  similar 
satisfactory  results,  see  Fig.  1  (right  column).  Since  this  represents  also  a 
globally-dependent  case,  we  have  performed  conventional  ICA  using  all  the 
pbcels,  not  to  our  surprise,  the  unexpected  overlaps  become  clearly  visible. 


Figure  4:  Normalized  pixel  TAG  projection  (left)  and  the  corresponding  factor 
TACs  (right),  obtained  by  VISDA  software. 

As  an  example  of  more  challenging  problems,  we  considered  dynamic  PET 
imaging  of  SERT  using  C^^McN"^5652.  Images  show  in  Fig.  5  (top  row) 
were  obtained  using  a  GE  4096  whole  body  scanner  with  15  slices  at  the  center 
'of  the  field-of-view.  The  slice  thickness  is  6.5  mm  and  the  spatial  resolution 
is  6''7  mm.  The  images  correspond  to  the  7th  PET  shce  and  were  taken 
between  8  to  18  time  snapshots  in  the  sequence.  Our  preliminary  results 
using  PICA,  shown  in  Fig.  5  (bottom  row),  were  quite  promising  in  that 
the  extracted  factors  closely  resemble  the  expected  compartmental  kinetics 
of  the  radioligand  consistent  with  its  pharmacology,  and  the  factor  images 
generated  simultaneously  reveal  regional  distribution  of  the  specific  (left)  and 
nonspecific  (right)  binding  sites.  We  believe  that  this  experiment  indicated 
a  relatively  successful  application  of  PICA  techmque  to  neuro- transporter 
binding  separation  in  PET,  given  the  difficulty  of  the  task. 

Final  experiment  reports  the  results  of  multivariate  clustering  in  blind 
estimation  of  both  the  input  function  and  factor  images  (microvessel  perme¬ 
ability  associated  with  fast  and  slow  perfusions-FP/SP).  The  separation  of 
plasma  space  and  fast  perfusion  permeability  is  expected  challenging.  We 
divided  the  TACs  into  two  portions  over  time.  We  firstly  estimated  FP  and 
SP  TACs  from  the  late  portion  of  the  TACs  assuming  insignificant  presence 
of  the  input  function.  We  then  removed  the  pixels  belonging  to  SP  and  es- 
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Figure  5:  Dynamic  PET  images  of  the  brain  and  the  reconstructed  factor  images 
that  represent  the  specific  and  non-specific  binding  sites  (2  biomarkers  of  interest). 


timated  the  input  function  from  the  early  portion  of  the  TACs  where  the 
FP  was  initialized  by  the  outcomes  of  the  preceding  step.  The  results  are 
given  in  Fig.  6.  The  processed  data  sets  are  insufficient  for  a  rigorous  test  of 
our  method.  So  far  as  we  can  tell,  it  is  roughly  compatible  with  the  clinical 
expectations,  but  it  must  be  regarded  as  unproved  until  it  has  been  checked 
against  more  realistic  truth. 
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FP  permeability 


SP  permeability 


Figure  6:  Blindly  estimated  input  function  &:  factor  images  by  pixel  TAG  clustering. 


DISCUSSION 


We  wish  to  suggest  that  our  preliminary  studies  indicated  a  promising  util¬ 
ity  of  hybrid  blind  source  separation  techniques  in  computed  simultaneous 


imaging  of  multiple  biomarkers.  Although,  the  optimality  of  this  method 
may  be  data-dependent,  we  would  expect  it  to  be  an  effective  tool  apphcable 
to  various  modalities.  We  are  currently  investigating  further  apphcation  of 
non-negative  matrix  factorization  [12]  and  separation  by  non-negative  reverse 
correlation  analysis,  where  the  centralldea  is  to  relax  bpth  the  independence 
and  nongaussian  assumptions. 
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ABSTRACT 

Dynamic  contrast  enhanced  magnetic  resonance  imaging  (DCE-MRI)  has  emerged  as  an  effective  tool  to  ac¬ 
cess  tumor  vascular  characteristics.  DCE-MRI  can  be  used  to  characterize  microvasculature  noninvasively  for 
providing  information  about  tumor  microvessel  structure  and  function  (e.g.,  tumor  blood  volume,  vascular  per¬ 
meability,  and  tumor  perfusion).  However,  pixels  of  DCE-MRI  represent  a  composite  of  more  than  one  distinct 
functional  biomarker  (e.g.,  microvessels  with  fast  or  slow  perfusion)  whose  spatial  distributions  are  often  het¬ 
erogeneous.  Complementary  to  various  existing  methods  (e.g.,  compartment  modeling,  factor  analysis),  this 
paper  proposes  a  blind  source  separation  method  that  allows  for  a  computed  simultaneous  imaging  of  multiple 
biomarkers  from  composite  DCE-MRI  sequences.  The  algorithm  is  based  on  a  partially-independent  component 
analysis,  whose  parameters  are  estimated  using  a  subset  of  informative  pixels  defining  the  independent  portion 
of  the  observations.  We  demonstrate  the  principle  of  the  approach  on  simulated  image  data  sets,  and  then  apply 
the  method  to  the  tissue  heterogeneity  characterization  of  breast  tumors.  As  a  result,  spatial  distribution  of 
tumor  blood  volume,  vascular  permeability,  and  tumor  perfusion,  as  well  as  their  time  activity  curves  (TACs) 
are  simultaneously  estimated. 

Keywords:  Independent  component  analysis  (ICA),  partially-independent  component  analysis  (PICA),  intrin¬ 
sic  dependency/non-intrinsic  dependency  of  the  components,  dynamic  contrast-enhanced  magnetic  resonance 
imaging  (DCE-MRI),  compartment  model,  time  activity  curves  (TACs). 

1.  INTRODUCTION 

Remarkable  advances  in  functional  imaging  have  been  made  in  developing  molecular-targeted  contrast  agents, 
ligands  and  imaging  probes.  Such  imaging  capabilities  will  allow  for  the  visualization  and  elucidation  of  im¬ 
portant  disease-causing  physiologic  and  molecular  processes  in  living  tissue.  Subsequently,  functional  imaging 
will  play  an  important  role  in  the  early  detection,  diagnosis,  and  treatment  of  diseases.^  It  is  known  that 
most  advanced  tumors  are  highly  heterogeneous  in  structure  that  may  reflect  the  underlying  angiogenesis  and/ or 
metastasis.^  Dynamic  contrast  enhanced  magnetic  resonance  imaging  (DCE-MRI)  is  a  noninvasive  imaging 
method  for  tumor  microvascular  characterization,  which  can  be  applied  to  assess  (and  potentially  predict)  the 
response  to  treatment  including  anti-angiogenic  drugs.  Kinetic  characteristics  changes  following  treatment  have 
correlated  with  histopathological  outcome  (e.g.,  microvessel  density)  and  patient  survival.  However,  widespread 
success  of  DCE-MRI  may  be  limited  by  the  need  for  further  technolog>^  development,  particularly  due  to  the 
lacking  of  quantitative  and  computational  data  analysis  tools  included  by  the  instruments. 

As  a  common  problem  in  functional  imaging,  pixels  represent  a  composite  of  more  than  one  distinct  molecular 
marker  (i.e.,  the  observed  pixel  intensity  will  consist  of  the  w'eighted  sum  of  activities  of  the  various  molecules). 
This  problem  exists  for  various  reasons,  e.g.,  target  mixture,  probe  non-specificity,  and  kinetics  or  spectrum 
overlap.  These  aspects  are  briefly  described  as  follows.  First,  mixed  signals  can  result  when  distinct  markers  are 
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combined  into  a  homogeneous  mixture  (e.g.,  fast  and  slow  flow  microvessels),  independent  of  spatial  resolution. 
Second,  ligand-receptor  binding  depends  largely  on  the  three-dimensional  shapes  of  both  elements,  where  a 
ligand  has  many  bonds  that  can  be  rotated  into  many  different  positions  resulting  in  many  shapes.  Third, 
even  with  a  precision  excitation  source,  any  overlap  of  the  absorption  spectra  of  the  fiuorophores  leads  to  the 
excitation  of  multiple  fiuorophores  whose  emission  spectra  often  also  overlap.  Thus,  the  observed  signal  intensity 
may  well  be  composed  of  the  emission  from  several  markers  of  differing  concentration  and  kinetics/spectrum 
(e.g.,  specific/nonspecific  bindings,  fast/slow  flows).  As  a  result,  the  overlap  of  multiple  molecular  signatures 
can  severely  decrease  the  sensitivity  and  specificity  for  the  measurement  of  molecular  signatures  associated  with 
difierent  disease  processes.  As  an  example,  imaging  neuro-transporters  in  the  brain  requires  the  passage  of 
radioligands  across  the  blood  brain  barrier  by  ways  of  their  high  lipophilicity.  But  hpophilicity  carries  the  risk 
of  high  nonspecific  binding  and  retention  in  the  white  matter  and  could  result  in  a  bias  of  the  estimated  kinetic 
parameters  that  are  used  to  measure  binding  to  specific  recognition  sites. 

It  is  well  known  that  Independent  Component  Analysis  (ICA)^^’^^  is  a  powerful  method  for  blind  source 
separation  with  a  strong  assumption  that  the  sources  are  independent  to  each  other.  This  paper  describes  a 
computation  approach  to  dependent  component  imaging,  where  functional  imaging  is  the  case.  The  method 
is  to  identify  an  informative  index  subspace  and  over  which  to  separate  mixed  imagery  sources  by  partially- 
independent  component  analysis  (PICA) ,  whose  parameters  are  estimated  using  informax  principle.  We  discuss 
the  theoretic  roadmap  of  the  approach,  and  its  applications  to  computer  simulation  phantoms  and  DCE-MRI 
sequences  of  breast  tumor. 


2.  THEORY  AND  METHOD 

Independent  component  analysis  (ICA)^®  is  a  statistical  and  computational  technique  for  revealing  hidden  factors 
that  underlie  sets  of  random  variables,  measurements,  or  signals.  The  application  of  ICA  has  been  found  in 
many  separate  fields  such  as  feature  extraction,  image  processing,  medical  image  processing,  telecommunication, 
econometric  signal  processing,  and  so  fourth. The  method  aims  at  recovering  the  unobservable  independent 
sources  (or  signals)  from  multiple  observed  data  masked  by  linear  or  nonhnear  mixing  of  the  components. One 
of  the  basic  assumptions  for  ICA  model  is  the  statistical  independence  between  components. However,  the 
dependent  components  are  often  occurred  in  the  real  world  situation,  including  functional  imaging  derived  from 
tissue  samples. 

2.1.  Compartment  Modeling 

Compartment  modeling  forms  the  basis  for  tracer  characterization  in  DCE-MRI.^  Fig.  1  shows  a  parallel  mode 
two-tissue  compartment  model.  ^  The  conventional  compartment  model  leads  to  a  set  of  first  order  differential 
equations: 

cy(t)  =  kifCp(t)-k2fCf(t) 

Cs(t)  =  kl^Cpit^  ^25^5  (t)  /-JN 

Ct{t)  =  Cf{t)  +  Cs{t)  ^ 

Cm  (t)  =  Cf  (t)  +  C5  (t)  +  Cp  (t) 

where  C/(t)  and  C5(t)  are  the  tissue  activity  in  the  fast  turnover  and  slow  turnover  pools,  respectively,  at  time 
t\  Cp(t)  is  the  tracer  concentration  in  plasma  (i.e.,  the  input  function);  Ct{t)  is  the  total  tissue  activity;  Cm(t) 
is  the  measured  total  tissue  activity;  kif  and  are  the  unidirectional  transport  constants  from  plasma  to 
tissue  (ml/min/g:  spatially  shift- varying);  and  k2f  and  k2$  are  the  rate  constants  for  effiux  (/min:  spatially 
shift- invariant).  It  is  important  to  note  that  Cj*(r),  C5(t),  Cp(t),  and  Cm{t)  are  also  called  the  time- activity  curves 
(TACs)  associated  with  a  pre-defined  region  of  interest  (ROI). 

It  can  be  shown  that  Cf  {t)  and  Cg  {t)  can  be  solved  analytically  in  a  parametric  form 

Cf{t)  =  ,2) 

Cs{t)  =  kisCp{t)(S^e 


where  (E>  denotes  the  mathematical  convolution  operation.  By  fitting  Cm(t)  to  the  measured  ROI  TAC  in  the 
light  of  pre- acquired  Cp(t)j  the  model  parameters  (A^i/,  kig^  /c2/,  k2s)  can  be  estimated.^’ ^ 
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Figure  1.  Two-tissue  compartmant  model  (parallel  mode). 


Based  on  linear  system  theory,  a  simple  method  can  be  developed  to  convert  temporal  kinetics  to  spatial 
information.^’^  First,  we  can  normalize  the  ROI  based  tissue  kinetics  to  define  three  TACs  for  each  pixel 
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(3) 


where  Vp  is  the  plasma  volume  in  tissue.  Second,  for  pixels  i  =  1, A?'  within  an  ROI,  we  let  kif[i)  and  kis{i) 
be  the  local  model  parameters  and  use  them  to  describe  the  dynamics  of  each  pixel  in  the  ROI 

=  kif{i)af{t)+kis{i)asit)  +Vp{i)ap{t)  (4) 

where  i)  is  the  measured  pixel  TAG,  kif(i)  and  kis{i)  are  the  permeability  of  fast  and  slow  turnover  regions 
in  the  pixel,  respectively,  and  Vp{i)  is  the  plasma  volume  in  the  pixel.  We  call  this  representation  as  factored 
compartment  modeling.^ 

Third,  let  t)e  the  sampling  time  points  of  the  DCE-MRI  measurements.  Then,  the  linear  least 

square  solution  of  Eq.  (4)  can  be  given  by  the  follow'ing  equation: 


(5) 


where 
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The  estimated  values  of  kif{i)^  kis{i)  and  Vp{i)  vary  from  pixel  to  pixel  and  reconstruct  three  factor  images 
respectively.  In  particular,  factor  images  kif{i)  and  kis{i)  represent  ROI  sub-regions  with  fast  and  slow  kinetics, 
respectively 

Preliminary  effort  has  been  recently  made  to  perform  blind  compartment  modeling  without  any  knowledge 
of  the  input  function. An  initial  effort,  eigenvector  based  multichannel  blind  deconvolution  (EVAM),^® 
was  used  to  estimate  the  parameters  of  a  two-tissue  compartment  model  for  PET  FDG  imaging, but  was 
shown  to  give  relatively  poor  (sensitive  to  noise)  and  non- unique  estimates  in  a  simulation  study. A  more 
optimal  solution  was  proposed  on  the  application  of  the  iterative  quadratic  maximum-likelihood  (IQML)  method 
to  parameter  estimation. The  blind  identification  problem  is  treated  as  a  nonlinear  least  square  problem 
whose  variables  are  separate.  Other  approaches  in  which  both  the  input  function  and  kinetic  parameters  are 
treated  as  unknowns  have  been  explored  in  [14,  15]. 


3,  PARTIALLY-INDEPENDENT  COMPONENT  ANALYSIS  (PICA) 

3.1.  Independent  Component  Analysis  (ICA) 

As  aforementioned,  one  potential  limitation  associated  with  compartment  analysis  is  that  they  are  all  restricted  to 
a  parametric  (thus  simplified)  model  that  may  not  adequately  describe  the  underlying  physiological  or  biochemical 
processes  about  tracer-target  interactions,  in  addition  to  the  likely  invasive  acquisition  of  the  input  function. 
Although  factor  analysis  (FA)  attempts  to  solve  the  problem,  the  results  were  mostly  unsatisfactory. 

From  hnear  system  theory, it  can  be  shown  that  the  solution  (zero-state  response)  to  a  kinetic  system  has 
the  very  general  form  as  shown  in  Eq.  (4),  or  E)q.  (7)  as  represented  in  vector -matrix  form.  This  motivates 
the  consideration  of  a  statistically-principled  computational  approach  involving  newly  invented  independent 
component  analysis  (ICA)  theory.^^’^®’^^  The  goal  is  to  blindly  and  computationally  reconstruct  both  A  and 
k  based  on  c^.  This  philosophy  for  computed  simultaneous  imaging  of  multiple  biomarkers  is  similar  in  spirit 
to  the  blind  source  separation  (BSS)  for  solving  the  cocktail-party  problem.^^ 

From  latent  variable  model  interpretation,^^  Eq.  (7) 

Cm(ht2) 


Cpi(ij  tn ) 

describes  how  the  observed  data  are  generated  by  a  process  of  mixing  the  latent  (or  “hidden”)  variables,  where 
matrix  A  is  called  the  mixing  matrix^  the  factor  images  (or  “source  signals”)  are  not  observable,  and  nothing  is 
known  about  the  properties  of  the  TACs  (or  “mixing  process”).  In  the  absence  of  this  information,  one  has  to 
proceed  “blindly”  to  recover  the  factor  images  from  their  TAG- modulated  activity  mixtures.^ 

We  can  state  such  computed  simultaneous  imaging  of  multiple  biomarkers  as  follows:  “Given  N  independent 
realizations  of  the  measured  pixel  TAG  vector  Cm(bt)i  i  find  an  estimate  of  the  inverse  of  the 

TAC-mixing  matrix  A(t)  and  factor  image  vector  k(^)  =  [A:i/(i), /ci5(^)]^-” 

ICA  method,  as  a  newly  invented  statistical  and  neural  computation  technique,  promises  a  powerful  computa¬ 
tional  tool  for  separating  hidden  sources  from  mixed  signals  when  many  classic  methods  fail  completely.^®  ICA 
method  utilizes  independence  as  a  guiding  principle  and  performs  BSS  based  on  a  nongaussian  factor  analysis 
with  a  unique  solution.^^  More  precisely,  by  assuming  that  the  hidden  components  are  statistically  independent 
with  nongaussian  distributions,  these  hidden  sources  can  be  found  by  ICA,  except  for  an  arbitrary  scaling  of 
each  signal  component  and  permutation  of  indices.  In  other  words,  it  is  feasible  to  find  a  demixing  matrix  W 
whose  individual  rows  are  a  rescaling  and  permutation  of  those  of  the  mixing  matrix  A.  ICA  approach  exploits 
primarily  temporal  diversity  in  that  the  dynamic  images  taken  at  different  times  carry  different  mixtures  of  the 
factor  images.®  There  are  several  algorithms  for  ICA  that  are  derived  from  different  optimization  principles. 
More  details  can  be  found  in  [11,  26]. 

3.2.  Partially-ICA 

We  have  found  that  direct  application  of  ICA  to  tumor  heterogeneity  characterization  using  all  the  pixels, 
however,  often  leads  to  an  unsatisfactory  recovery  of  factor  images  k(i).  By  a  closer  look  at  the  joint  distribution 
of  the  factor  images,  we  found  that  they  are  often  not  statistically  independent  over  the  whole  pixel  set.®  This 
shall  not  be  a  surprise  since  factor  images  are  expected  to  be  piece- wise  continuous  thus  form  clusters  over  the 
joint  distribution.  It  can  be  further  concluded  that  such  joint  distribution  clusters  correspond  to  the  overlapped 
homogeneous  areas  of  the  factor  images.  Thus,  we  shall  expect  to  achieve  a  better  factor  image  decomposition 
using  a  subset  of  pixels  that  supports  the  independency  of  the  factor  images. 

Inspired  by  such  reasoning,  we  proposed  a  partially-ICA  (PICA)  technique  in  [5].  Rather  than  using  all 
the  pixels  that  give  rise  to  a  large  decomposition  error  due  to  source  dependency,  we  attempt  to  (iteratively) 


=  A  kuii)  ,  (7) 


identify  a  pixel  subset  supporting  source  independency  and  over  which  to  estimate  the  demixing  matrix  W*  and 
subsequently  factor  images  k(i). 

Compared  with  the  basic  ICA  model,  where  each  observation  is  a  linear  combination  of  independent  compo¬ 
nents,  our  PICA  model  assumes  that  each  observation  Xi  is  a  linear  mixture  of  statistically  dependent  components 
with  an  n  by  n  non-singular  mixing  matrix  A,  i.e., 
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and  m  is  the  number  of  pixels  in  each  functional  image.  Our  task  is  to  recover  the  dependent  components  from 
observations  by  still  utilizing  ICA. 


Let  us  review  the  mechanism  of  independent  component  analysis  on  basic  ICA  model.  The  components  Si 
are  statistically  independent,  while  based  on  the  Central  Limit  Theorem;  the  distribution  of  a  sum  (observations) 
of  independent  random  variables  (components)  tends  toward  a  Gaussian  distribution,  under  certain  conditions. 
Therefore,  the  procedure  that  ICA  searches  for  estimates  of  the  components  is  to  find  directions,  such  that 
the  projections  of  observations  on  each  direction  are  distributed  with  most  non- Gaussian  distribution.  This 
is  the  reason  that  ICA  algorithm  can  help  find  the  (statistically)  precise  estimate  of  the  components  if  the 
components  are  completely  independent  (except  those  two  ambiguities  of  ICA  on  the  scale  and  order  of  the 
components) .  However,  it  will  mislead  direction  finding  if  the  components  are  dependent  but  ICA  algorithm  is 
still  superimposed  on  the  observations,  which  are  mixtures  of  the  dependent  components.  Also,  it  will  give  rise 
to  a  large  separation  error  since  all  the  pixels  are  utilized  for  ICA  calculation  while  the  components  over  all  these 
pixels  are  statistically  dependent. 

In  order  to  effectively  utilize  ICA  for  this  problem,  an  informative  pixel  index  subspace  (corresponding  to  the 
independent  part  of  the  components)  needs  to  be  identified.  Then  we  can  perform  ICA  over  this  subspace  to 
recover  the  estimated  mixing  matrix  over  the  subspace.  Imposing  this  mixing  matrix  over  all  the  pixel  indices, 
the  dependent  components  are  then  available  to  be  recovered.  The  key  point  is  to  identify  the  informative  pixel 
indices,  over  which  the  components  are  statistically  independent.  The  difficulty  of  this  approach  is  that  the 
independent  subspace  of  the  components  needs  to  be  identified  without  any  statistical  information  derived  from 
the  components  themselves:  the  components  are  just  what  need  to  be  recovered.  We  only  have  information 
derived  from  observations  rather  than  from  components.  For  simplicity,  we  consider  the  dependency  of  the 
components  rather  than  that  of  the  observations.  It  is  well  know  that  the  statistical  dependency/independency 
of  the  components  s  could  be  measured  (visualized)  by  means  of  scatter  plot  with  m  sample  points  iS'i,  ^2, 
in  n  dimensional  space  (each  dimension  corresponds  to  a  component):  they  consist  of  dependent/independent 
components,  which  we  can  perform  linear/nonlinear  regression  curve  to  estimate  the  statistical  relation  between 
n  components. 

Clearly,  when  the  number  of  pixels  m  is  much  larger  than  the  number  of  components  n,  w'hich  is  the  situation 
for  functional  imaging,  most  probably  the  sample  points  are  clustered  into  some  clusters.  Based  on  this  obser¬ 
vation,  we  divide  the  dependency  of  the  components  into  two  categories:  intrinsic  dependency  and  non-intrinsic 
dependency.  By  intrinsic  dependency  of  the  components,  we  mean  the  dependency  caused  by  the  linear  and/or 
nonlinear  correlation  between  components  over  cluster  centers.  We  refer  to  non-intrinsic  dependency  of  the 
components  as  the  dependency  over  sample  points  inside  each  cluster.  In  other  words,  the  intrinsic  dependency 
corresponds  to  the  global  dependency  among  clusters,  while  non-intrinsic  dependency  corresponds  to  the  local 


dependency  among  samples  in  each  cluster.  For  ICA  to  work  well,  we  need  to  remove  both  intrinsic  dependency 
and  non-intrinsic  dependency  of  the  components.  This  can  be  accomplished  by  removing  the  pixels  that  con¬ 
tribute  to  intrinsic  dependency  and  non-intrinsic  dependency  of  the  components,  while  retaining  the  informative 
pixel  indices.  Finally,  we  will  apply  the  ICA  onto  the  subset  of  pixels  for  effective  recovery  of  the  components 
over  those  informative  pixel  indices. 

Notice  that  both  intrinsic  dependency  and  non-intrinsic  dependency  of  the  components  should  be  removed 
with  only  the  information  from  observations  rather  than  from  components.  A  possible  approach  to  remove  intrin¬ 
sic  and  non-intrinsic  dependencies  can  be  summarized  as  follows.  The  intrinsic  dependency  of  the  components 
can  be  removed  by  removing  all  of  the  clusters  except  one.  The  non-intrinsic  dependency  of  the  components  can 
be  alleviated  by  removing  some  sample  points  in  the  remained  cluster.  We  will  describe  the  removal  procedures 
next. 

3.2.1.  Intrinsic  dependency  removal 

In  order  to  remove  intrinsic  dependency  of  the  components  from  the  observations,  the  joint  distribution  of  the 
observations  is  estimated  with  Expectation-Maximization  (EM)  algorithm  initialized  by  k-means  method.  We 
assume  that  the  number  of  clusters  t  is  known  with  some  prior  knowledge  of  the  problem.  In  our  DCE-MRI 
study  of  breast  tumor,  the  prior  knowledge  tells  us  that  this  number  should  be  two  considering  the  fast-flow  and 
slow-flow  characteristics  of  the  breast  tumor.  Assume  that  the  distribution  of  the  observations  is  in  the  following 
form: 


t 

where  p{x;  is  the  ith  Gaussian  distribution  with  and  af  as  its  mean  and  variance,  and  tti  is  the  weight 

t 

of  the  zth  Gaussian  distribution,  =  1.  Notice  that  each  Gaussian  distribution  forms  a  cluster  in  scatter 

i=:l 

plot  of  the  observations. 

We  remove  intrinsic  dependency  of  the  components  in  a  statistical  way.  Each  sample  point  Xj  is  removed  in 
observation  scatter  plot  statistically,  with  the  probability  of 

of  X,) = (11) 

where  k  is  the  index  of  the  remaining  cluster. 

3.2.2.  Non-intrinsic  dependency  Removal 

For  removal  of  non-intrinsic  dependency  of  the  components,  the  remained  sample  points  in  scatter  plot  of  the 
observations  are  statistically  down-sampled  with  Pe  as  its  parameter,  i.e.,  Xj  is  removed  statistically  with  the 
probability  of 


As  a  result,  the  remained  sample  points  in  scatter  plot  of  the  observations  w^ould  be  uniformly  distributed.  It  is 
important  that  the  statistical  dependency  of  the  components  over  the  corresponding  sample  points  are  removed 
without  any  destruction  of  the  linear  dependency  of  the  observations  over  the  remained  sample  points.  In  fact, 
this  dependency  should  not  be  removed  for  the  following  ICA  computation,  because  the  observations  over  the 
remained  sample  points  are  still  the  linear  mixing  of  the  components,  which  are  statistically  independent  over 
the  corresponding  sample  points  from  the  PICA  model. 


Observation  Intrinsic  dependency  removal  Remaining  Cluster' 


Figure  2.  Removal  of  the  intrinsic  dependency  between  the  components:  (i)  observation  consists  of  more  than  one  cluster; 
(ii)  the  removal  of  intrinsic  cluster  dependency;  (iii)  the  remaining  cluster  after  the  removal  of  intrinsic  dependency. 


3.2.3.  Proposed  Algorithm 

Assume  I  is  the  remained  sample  point  index  subset,  which  is  obtained  by  the  removal  of  intrinsic  dependency 
and  non-intrinsic  dependency  of  the  components  from  the  scatter  plot  of  the  observations.  I  corresponds  to  the 
independent  pixel  index  subspace  of  the  components.  We  have  following  procedure  for  a  composite  separation 
of  observations: 

Step  1:  estimate  the  joint  distribution  p{x)  of  the  observations  with  EM  algorithm  initialized  by  k-means 
method,  where  the  number  of  clusters  t  is  assumed  to  be  known  from  some  prior  knowledge; 

Step  2:  remove  intrinsic  dependency  of  the  components  by  a  statistical  cluster  removal  method  with  eq.  11, 
and  remove  non-intrinsic  dependency  of  the  components  by  a  statistical  down-sampling  method  by  utilizing  eq. 
12  with  parameter  P^,  both  in  the  scatter  plot  of  the  observations;  assume  the  retained  pixel  index  subset  is  /; 

Step  3:  perform  ICA  on  the  retained  part  of  the  observations,  {Xj^j  €  /},  to  obtain  the  estimated  mixing 
matrix  A.  Since  the  observations  over  the  retained  pixel  index  subspace  I  are  independent,  which  satisfies  the 
assumption  of  basic  ICA  model;  the  mixing  matrix  A  is  expected  to  be  better  estimated  except  the  ambiguity 
of  scaling  and  ordering. 

Step  4:  impose  the  estimated  mixing  matrix  on  the  observations  over  entire  indices  to  obtain  the  estimated 
components,  i.e.,  s  =  A”^x. 

There  are  two  parameters,  t  and  in  the  above  algorithm,  t  is  set  with  some  number  by  prior  knowledge 
of  the  problem,  while  Pe  is  a  trade-off  parameter  for  controlling  the  uniformity  of  the  distribution  over  the 
sparseness  of  the  sample  points  in  the  scatter  plot  of  observations  after  the  removal  procedures.  The  smaller  the 
Pe  is,  the  more  uniform  the  distribution  of  sample  points  over  the  index  subspace  looks  like,  and  the  less  sample 
points  over  the  index  subspace  will  be.  The  directions  searched  out  for  the  sample  points  over  the  index  subspace 
with  ICA  algorithm  determine  the  rows  of  the  estimated  mixing  matrix.  Theoretically  speaking,  if  the  parameter 
Pe  is  set  smaller,  the  more  uniformed  distribution  of  the  sample  points  in  scatter  plot  of  the  observations  over 
the  index  subspace  will  merit  this  direction  searching;  On  the  other  hand,  the  number  of  sample  points  becomes 
less  that  limits  the  sample  points  to  form. that  distribution,  demerging  this  direction  searching. 

4.  RESULT  AND  DISCUSSION 

We  first  applied  our  method  to  two  computer  simulation  phantoms  for  the  removal  of  intrinsic  dependency 
and  non-intrinsic  dependency  among  components.  Then  we  applied  our  method  to  a  data  set  generated  by- 
compartment  models,  and  a  real  DCE-MRI  data  set  of  tumor  heterogeneity  study. 


(b)  The  situation  of  non-linear  dependency  between  components 


Figure  3.  The  scatter  plots  of  the  component  recovery  process,  for  both  linear  (a)  and  non-linear  (b)  dependency 
cases,  with  ICA  and  PICA  method,  (i)  components,  (ii)  observations,  (iii)  recovered  components  from  ICA  method,  (iv) 
observations  over  the  informative  index  subspace,  (v)  recovered  components  over  the  informative  index  subspace,  (vi) 
recovered  components  over  full  index  space  from  PICA  method. 

4.1.  Phantom  study 

Figure  2(i)  shows  the  data  set  that  is  generated  by  a  sum  of  two  Gaussian  distributions:'  centered  at  (0,1)  and 
(1,0)  with  standard  variations  as  two  components.  Clearl}^  the  components  are  not  independent  (i.e.,  they  are 
intrinsically  dependent),  and  the  removal  of  either  Gaussian  distribution  makes  the  sample  points  over  the  index 
subspace  statistically  independent  (see  Figure  2(ii)  and  (iii)). 

The  next  data  set  is  generated  such  that  the  two  components  are  independent  in  their  first  half  (left  half), 
and  have  linear/non-linear  correlation  for  the  other  half  (right  half),  both  with  uniformly  distributed  intensities. 
The  mixing  matrix  is  randomly  generated  to  form  observations.  Clearly  there  is  a  linear/non-linear  non-intrinsic 
dependency  between  the  components  in  this  experiment.  Figure  3  (upper/lower  figure)  shows  the  scatter  plots  in 
the  component  recovery  process  with  ICA  and  PICA  methods  for  the  situation  of  linear/non-linear  dependency 
between  the  components.  Figure  4  shows  the  component  recovery  process,  where  the  randomized  intensities  of  the 
first  half  and  the  second  half  are  shown  in  two-dimensional  images.  Note  that  the  pixels  for  the  first  component 
in  the  figure  are  reordered  according  to  its  intensities.  Evidently,  the  estimated  components  recovered  by  PICA 
method  are  much  closer  to  the  ground  truth  of  the  components  by  comparing  both  the  scatter  plots  in  Figure 
3(i),  (iii),  (vi)  and  the  images  in  Figure  4(ii),  (iv),  (vi). 

4.2.  Experiments  on  compartment  model 

A  compartment  model  is  used  to  simulate  the  dynamic  behavior  of  the  breast  tumor  obtained  with  a  DCE- 
MRI  sequence.  The  mask  for  the  fast-flow  (FF) /slow-flow  (SF)  patterns,  as  well  as  the  overlap  region  of  FF 
and  SF,  is  shown  in  Figure  5(a).  The  pixel  intensity  in  the  FF-dominant  region/overlap  region/SF-dominant 
region  is  a  linear  combination  of  FF  and  SF  patterns  with  the  corr^ponding  intensity  weights  of  0.9/0.5/0.1 
and  0.1/0.5/0.9.  In  Figure  5(a),  bright/dark  gray  color  corresponds  to  FF/SF  region,  and  the  light  gray/dark 
corresponds  to  overlap/ background,  respectively.  FVom  the  compartment  model,  w'e  can  get  a  sequence  of  images 
shown  in  Figure  5(b).  Our  task  is  to  identify  the  FF  and  SF  patterns  hidden  in  the  observed  sequence  of  images. 


(i) 


(Vi) 


(ii)  (iii)  (iv)  (v) 

(a)  The  situation  of  linear  dependency  between  components 


(i)  (ii)  (iii)  (iv)  (v)  (vi) 

(b)  The  situation  of  non-linear  dependency  between  components 


Figure  4.  Recovery  results  from  ICA  and  PICA,  (i)  observations,  (ii)  recovered  components  by  utilizing  ICA  method, 
(iii)  observations  over  the  informative  index  subspace,  (iv)  recovered  components  over  the  informative  index  subspace, 
(v)  recovered  components  over  full  index  space  by  utilizing  PICA  method,  (vi)  real  components. 


(a)  (b) 


Figure  5 .  A  simulater  tumor  phantom  inclusing  fast  and  slow  kinetic  subregions;  (a)  Image  mzisk  in  compartment  model, 
(b)  The  sequence  of  images  from  the  compartment  model  with  the  mask  shown  in  (a) 


We  applied  our  PICA  method  to  separate  the  observed  composite  images  into  two  FF  and  SF  patterns,  as  well 
as  to  have  the  time-activity  curve  (TAG).  The  only  two  spatial  patterns,  FF  and  SF,  which  is  a  prior  knowledge 
of  this  problem,  induced  us  to  make  an  assumption  of  two  components  from  the  observed  image  sequence.  Thus, 
we  divided  the  image  sequence  into  the  first  half  of  the  sequence  and  the  second  half  of  the  sequence,  according 
to  the  time  index.  Then  PICA  was  performed  for  each  pair  of  observed  images,  one  from  the  first  half,  and 
another  from  the  second  half,  and  the  corresponding  components  were  estimated.  The  final  estimation  of  the 
components  (spatial  FF/SF  patterns)  was  attained  by  averaging  all  of  the  estimated  components  obtained  from 
each  pair  of  observed  images.  Figure  6  shows  the  spatial  FF/SF  pattern  and  the  corresponding  TACs  achieved 
by  PICA,  together  wdth  those  by  ICA  for  comparison.  Notice  that  each  TAG  according  to  the  compartment 
model  should  be  a  positive  curve.  Figure  6(a)  shows  a  better  performance  for  the  estimation  of  the  FF  and 


SF  patterns,  which  are  closer  to  the  ground  truth  of  0.9;0.5:0.1  and  0.1:0.5:0.9/0.9:0.3:0.1  and  0.1:0.7;0.9  from 
Figure  6(a)(i)/(ii).  It  is  shown  from  Figure  6(b)  that  TACs  from  PICA  are  better  than  those  from  ICA,  since 
they  are  both  approximately  positive  and  more  fit  to  the  meaning  of  dynamic  fast-fiow  and  slow-flow  activity  of 
the  tumor  tissue. 


ICA  PICA  ICA  PICA 


(a)  The  recovered  spatial  FF/SF  patterns 


(b)  The  corresponding  TACs 


Figure  6.  ICA  and  PICA  result,  (a)  the  factor  images  results  for  compartment  model  with  FF:overlap:SF  region  inten¬ 
sity  for  FF  pattern  and  for  SF  pattern  to  be  0.9:0.5:0.1  and  0.1:0.5:0.9/0.9:0.3:0.1  and  0.1:0.7:0.9  respectively;  (b)  the 
corrsepondiag  TACs. 


4.3.  Experiment  on  real  DCE-MRI  data  set 

We  tested  our  PICA  method  with  a  real  DCE-MRI  sequence  of  breast  tumor  studies.  Figure  7  shows  a  typical 
sequence  of  breast  tumor  DCE-MRI  study.  Compared  with  results  from  the  direct  application  of  ICA,  our  results 
using  PICA  shown  in  Figure  8,  were  quite  promising  in  the  extracted  factors  that  closely  resemble  the  expected 
characteristics  of  compartmental  kinetics  of  tumors.  The  factor  images  and  the  TACs  reveal  regional  distribution 
of  the  FF  and  SF  patterns. 
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Figure  7.  (a)  Tumor  (breast  cancer)  spatial  heterogeneity  revealed  by  DCE-MRI;  (b)  Sequences  of  the  breast  tumor 
ROl 

ICA  PICA  TACs  from  PICA 


Figure  8.  The  recovered  factor  images  of  the  breast  cancer  spatial  heterogeneity  and  the  corresponding  TACs  revealed 
by  DCE-MRI  with  PICA  method 
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